PR data

Import data

data <- read.csv("output/pr_curve_extracted_rates.csv")
str(data)
## 'data.frame':    1114 obs. of  9 variables:
##  $ X                : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Species          : chr  "P. grandis" "P. grandis" "P. grandis" "P. grandis" ...
##  $ colony_id.x      : chr  "POG-150" "POG-150" "POG-150" "POG-150" ...
##  $ avg_temp_interval: num  17.2 17.1 21.1 21.1 24.1 ...
##  $ Temperature      : int  17 17 21 21 24 24 28 28 30 30 ...
##  $ Light            : int  0 560 0 560 0 560 0 560 0 560 ...
##  $ Run              : int  4 4 4 4 4 4 4 4 4 4 ...
##  $ micromol.cm2.s   : num  -0.0583 0.0334 -0.0762 0.0421 -0.0981 ...
##  $ micromol.cm2.h   : num  -210 120 -274 152 -353 ...
data$Temperature <- as.numeric(data$avg_temp_interval)
str(data)
## 'data.frame':    1114 obs. of  9 variables:
##  $ X                : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Species          : chr  "P. grandis" "P. grandis" "P. grandis" "P. grandis" ...
##  $ colony_id.x      : chr  "POG-150" "POG-150" "POG-150" "POG-150" ...
##  $ avg_temp_interval: num  17.2 17.1 21.1 21.1 24.1 ...
##  $ Temperature      : num  17.2 17.1 21.1 21.1 24.1 ...
##  $ Light            : int  0 560 0 560 0 560 0 560 0 560 ...
##  $ Run              : int  4 4 4 4 4 4 4 4 4 4 ...
##  $ micromol.cm2.s   : num  -0.0583 0.0334 -0.0762 0.0421 -0.0981 ...
##  $ micromol.cm2.h   : num  -210 120 -274 152 -353 ...
data <- data %>%
  arrange(colony_id.x) %>%
  group_by(colony_id.x) %>%
  mutate(curve_id = group_indices())


data %>% 
ggplot(aes(x=Temperature, y=micromol.cm2.s, color=Species, group=colony_id.x))+ 
  geom_point()+
  geom_line()+
  ylab("µmol O2 cm^2 s-1")+
  facet_wrap(~ Species*Light,  ncol = 5)+
  theme(legend.position = "none")

d <- data %>%
  filter(Light==0) %>%
  filter(!Species == "Recruit")

d$micromol.cm2.s <- -d$micromol.cm2.s
d$micromol.cm2.s <- replace(d$micromol.cm2.s, d$micromol.cm2.s<0,0)

d <- d[,-c(1,4,6,7,9)]
str(d)
## gropd_df [278 × 5] (S3: grouped_df/tbl_df/tbl/data.frame)
##  $ Species       : chr [1:278] "P. grandis" "P. grandis" "P. grandis" "P. grandis" ...
##  $ colony_id.x   : chr [1:278] "POG-150" "POG-150" "POG-150" "POG-150" ...
##  $ Temperature   : num [1:278] 17.2 21.1 24.1 28 29.7 ...
##  $ micromol.cm2.s: num [1:278] 0.0583 0.0762 0.0981 0.1325 0.1434 ...
##  $ curve_id      : int [1:278] 1 1 1 1 1 1 1 1 1 2 ...
##  - attr(*, "groups")= tibble [30 × 2] (S3: tbl_df/tbl/data.frame)
##   ..$ colony_id.x: chr [1:30] "POG-150" "POG-151" "POG-156" "POG-159" ...
##   ..$ .rows      : list<int> [1:30] 
##   .. ..$ : int [1:9] 1 2 3 4 5 6 7 8 9
##   .. ..$ : int [1:9] 10 11 12 13 14 15 16 17 18
##   .. ..$ : int [1:8] 19 20 21 22 23 24 25 26
##   .. ..$ : int [1:8] 27 28 29 30 31 32 33 34
##   .. ..$ : int [1:9] 35 36 37 38 39 40 41 42 43
##   .. ..$ : int [1:8] 44 45 46 47 48 49 50 51
##   .. ..$ : int [1:9] 52 53 54 55 56 57 58 59 60
##   .. ..$ : int [1:9] 61 62 63 64 65 66 67 68 69
##   .. ..$ : int [1:9] 70 71 72 73 74 75 76 77 78
##   .. ..$ : int [1:8] 79 80 81 82 83 84 85 86
##   .. ..$ : int [1:8] 87 88 89 90 91 92 93 94
##   .. ..$ : int [1:8] 95 96 97 98 99 100 101 102
##   .. ..$ : int [1:8] 103 104 105 106 107 108 109 110
##   .. ..$ : int [1:9] 111 112 113 114 115 116 117 118 119
##   .. ..$ : int [1:8] 120 121 122 123 124 125 126 127
##   .. ..$ : int [1:9] 128 129 130 131 132 133 134 135 136
##   .. ..$ : int [1:9] 137 138 139 140 141 142 143 144 145
##   .. ..$ : int [1:8] 146 147 148 149 150 151 152 153
##   .. ..$ : int [1:8] 154 155 156 157 158 159 160 161
##   .. ..$ : int [1:9] 162 163 164 165 166 167 168 169 170
##   .. ..$ : int [1:9] 171 172 173 174 175 176 177 178 179
##   .. ..$ : int [1:9] 180 181 182 183 184 185 186 187 188
##   .. ..$ : int [1:7] 189 190 191 192 193 194 195
##   .. ..$ : int [1:7] 196 197 198 199 200 201 202
##   .. ..$ : int [1:32] 203 204 205 206 207 208 209 210 211 212 ...
##   .. ..$ : int [1:9] 235 236 237 238 239 240 241 242 243
##   .. ..$ : int [1:9] 244 245 246 247 248 249 250 251 252
##   .. ..$ : int [1:9] 253 254 255 256 257 258 259 260 261
##   .. ..$ : int [1:8] 262 263 264 265 266 267 268 269
##   .. ..$ : int [1:9] 270 271 272 273 274 275 276 277 278
##   .. ..@ ptype: int(0) 
##   ..- attr(*, ".drop")= logi TRUE
colnames(d) <- c("Species", "colony_id", "temp","rate","curve_id")
d$curve_id <- as.numeric(d$curve_id)
str(d)
## gropd_df [278 × 5] (S3: grouped_df/tbl_df/tbl/data.frame)
##  $ Species  : chr [1:278] "P. grandis" "P. grandis" "P. grandis" "P. grandis" ...
##  $ colony_id: chr [1:278] "POG-150" "POG-150" "POG-150" "POG-150" ...
##  $ temp     : num [1:278] 17.2 21.1 24.1 28 29.7 ...
##  $ rate     : num [1:278] 0.0583 0.0762 0.0981 0.1325 0.1434 ...
##  $ curve_id : num [1:278] 1 1 1 1 1 1 1 1 1 2 ...
##  - attr(*, "groups")= tibble [30 × 2] (S3: tbl_df/tbl/data.frame)
##   ..$ colony_id: chr [1:30] "POG-150" "POG-151" "POG-156" "POG-159" ...
##   ..$ .rows    : list<int> [1:30] 
##   .. ..$ : int [1:9] 1 2 3 4 5 6 7 8 9
##   .. ..$ : int [1:9] 10 11 12 13 14 15 16 17 18
##   .. ..$ : int [1:8] 19 20 21 22 23 24 25 26
##   .. ..$ : int [1:8] 27 28 29 30 31 32 33 34
##   .. ..$ : int [1:9] 35 36 37 38 39 40 41 42 43
##   .. ..$ : int [1:8] 44 45 46 47 48 49 50 51
##   .. ..$ : int [1:9] 52 53 54 55 56 57 58 59 60
##   .. ..$ : int [1:9] 61 62 63 64 65 66 67 68 69
##   .. ..$ : int [1:9] 70 71 72 73 74 75 76 77 78
##   .. ..$ : int [1:8] 79 80 81 82 83 84 85 86
##   .. ..$ : int [1:8] 87 88 89 90 91 92 93 94
##   .. ..$ : int [1:8] 95 96 97 98 99 100 101 102
##   .. ..$ : int [1:8] 103 104 105 106 107 108 109 110
##   .. ..$ : int [1:9] 111 112 113 114 115 116 117 118 119
##   .. ..$ : int [1:8] 120 121 122 123 124 125 126 127
##   .. ..$ : int [1:9] 128 129 130 131 132 133 134 135 136
##   .. ..$ : int [1:9] 137 138 139 140 141 142 143 144 145
##   .. ..$ : int [1:8] 146 147 148 149 150 151 152 153
##   .. ..$ : int [1:8] 154 155 156 157 158 159 160 161
##   .. ..$ : int [1:9] 162 163 164 165 166 167 168 169 170
##   .. ..$ : int [1:9] 171 172 173 174 175 176 177 178 179
##   .. ..$ : int [1:9] 180 181 182 183 184 185 186 187 188
##   .. ..$ : int [1:7] 189 190 191 192 193 194 195
##   .. ..$ : int [1:7] 196 197 198 199 200 201 202
##   .. ..$ : int [1:32] 203 204 205 206 207 208 209 210 211 212 ...
##   .. ..$ : int [1:9] 235 236 237 238 239 240 241 242 243
##   .. ..$ : int [1:9] 244 245 246 247 248 249 250 251 252
##   .. ..$ : int [1:9] 253 254 255 256 257 258 259 260 261
##   .. ..$ : int [1:8] 262 263 264 265 266 267 268 269
##   .. ..$ : int [1:9] 270 271 272 273 274 275 276 277 278
##   .. ..@ ptype: int(0) 
##   ..- attr(*, ".drop")= logi TRUE

TPC fitting Padifeld et al rTPC and nls.multstart: A new pipeline to fit thermal performance curves in r
https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13585

Sharpe Schoolfield 1981 model Schoolfield, R. M., Sharpe, P. J. H., & Magnuson, C. E. (1981). Non-linear regression of biological temperature-dependent rate models based on absolute reaction-rate theory. Journal of theoretical biology, 88(4), 719-731. https://doi.org/10.1016/0022-5193(81)90246-0

#edit nls_multstart to allow for a progress bar

# edit nls_multstart to allow for a progress bar
nls_multstart_progress <- function(formula, data = parent.frame(), iter, start_lower, 
                                   start_upper, supp_errors = c("Y", "N"), convergence_count = 100, 
                                   control, modelweights, ...){
  if(!is.null(pb)){
    pb$tick()
  }
  nls_multstart(formula = formula, data = data, iter = iter, start_lower = start_lower, 
                start_upper = start_upper, supp_errors = supp_errors, convergence_count = convergence_count, 
                control = control, modelweights = modelweights, ...)
}

start progress bar and estimate time it will take

# start progress bar and estimate time it will take
number_of_models <- 1
number_of_curves <- length(unique(d$curve_id))

# setup progress bar
pb <- progress::progress_bar$new(total = number_of_curves*number_of_models,
                                 clear = FALSE,
                                 format ="[:bar] :percent :elapsedfull")

fit chosen model formulation in rTPC

# fit chosen model formulation in rTPC
d_fits <- nest(d, data = c(temp, rate)) %>%
  mutate(sharpeschoolhigh = map(data, ~nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = .x,
                        iter = c(3,3,3,3),
                        start_lower = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981') - 1,
                        start_upper = get_start_vals(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981') + 1,
                        lower = get_lower_lims(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981'),
                        upper = get_upper_lims(.x$temp, .x$rate, model_name = 'sharpeschoolhigh_1981'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)))
d_preds <- mutate(d_fits, new_data = map(data, ~tibble(temp = seq(min(.x$temp), max(.x$temp), length.out = 100)))) %>%
  # get rid of original data column
  select(., -data) %>%
  # stack models into a single column, with an id column for model_name
  pivot_longer(., names_to = 'model_name', values_to = 'fit', c(sharpeschoolhigh)) %>%
  # create new list column containing the predictions
  # this uses both fit and new_data list columns
  mutate(preds = map2(fit, new_data, ~augment(.x, newdata = .y))) %>%
  # select only the columns we want to keep
  select(curve_id, Species,model_name, preds) %>%
  # unlist the preds list column
  unnest(preds)

glimpse(d_preds)
## Rows: 3,000
## Columns: 6
## Groups: colony_id [30]
## $ colony_id  <chr> "POG-150", "POG-150", "POG-150", "POG-150", "POG-150", "POG…
## $ curve_id   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ Species    <chr> "P. grandis", "P. grandis", "P. grandis", "P. grandis", "P.…
## $ model_name <chr> "sharpeschoolhigh", "sharpeschoolhigh", "sharpeschoolhigh",…
## $ temp       <dbl> 17.21817, 17.44516, 17.67214, 17.89913, 18.12611, 18.35309,…
## $ .fitted    <dbl> 0.05791013, 0.05903478, 0.06017729, 0.06133772, 0.06251611,…
ggplot(d_preds) +
  geom_line(aes(temp, .fitted, col = Species, group=curve_id)) +
  geom_point(aes(temp, rate, col = Species), d) +
  #facet_wrap(~curve_id*Species, scales = 'free_y', ncol = 6) +
  theme_bw() +
  #theme(legend.position = 'none') +
  scale_color_brewer(type = 'qual', palette = 2) +
  labs(x = 'Temperature (ºC)',
       y = 'Light Enhanced Dark Respiration rate',
       title = 'sharpeschoolhigh')

ggplot(d_preds) +
  geom_line(aes(temp, .fitted, col = Species, group=curve_id)) +
  geom_point(aes(temp, rate, col = Species), d) +
  facet_wrap(~Species,  ncol = 4) +
  theme_bw() +
  theme(legend.position = 'none') +
  scale_color_brewer(type = 'qual', palette = 2) +
  labs(x = 'Temperature (ºC)',
       y = 'Light Enhanced Dark Respiration rate',
       title = 'sharpeschoolhigh')

#calculate TPC metrics by species

d_params <- pivot_longer(d_fits, names_to = 'model_name', values_to = 'fit', c(sharpeschoolhigh)) %>%
  mutate(params = map(fit, calc_params)) %>%
  select(curve_id, Species, model_name, params) %>%
  unnest(params)

glimpse(d_params)
## Rows: 30
## Columns: 15
## Groups: colony_id [30]
## $ colony_id             <chr> "POG-150", "POG-151", "POG-156", "POG-159", "POG…
## $ curve_id              <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1…
## $ Species               <chr> "P. grandis", "P. grandis", "P. grandis", "P. gr…
## $ model_name            <chr> "sharpeschoolhigh", "sharpeschoolhigh", "sharpes…
## $ rmax                  <dbl> 0.1341227, 0.1692973, 0.1834050, 0.1645734, 0.14…
## $ topt                  <dbl> 31.46017, 32.99975, 33.07686, 32.19303, 33.44543…
## $ ctmin                 <dbl> -5.66082708, -5.93625000, -5.69314000, -9.320970…
## $ ctmax                 <dbl> 55.55617, 48.59375, 44.51986, 54.66203, 47.32543…
## $ e                     <dbl> 0.5622984, 0.5017270, 0.4942046, 0.3974826, 0.48…
## $ eh                    <dbl> 0.3230061, 0.9464816, 1.6929276, NA, 0.9786143, …
## $ q10                   <dbl> 2.088531, 1.917583, 1.890915, 1.675354, 1.883673…
## $ thermal_safety_margin <dbl> 24.096, 15.594, 11.443, 22.469, 13.880, 10.085, …
## $ thermal_tolerance     <dbl> 61.217, 54.530, 50.213, 63.983, 54.726, 48.264, …
## $ breadth               <dbl> 11.432, 9.517, 8.169, 11.732, 9.195, 7.654, 8.30…
## $ skewness              <dbl> 0.23929229, -0.44475468, -1.19872297, NA, -0.490…
TPC.pars <- as.data.frame(d_params)
summary(aov(rmax ~ Species, TPC.pars))
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## Species      3 0.001081 0.0003603   0.868   0.47
## Residuals   26 0.010789 0.0004150
summary(aov(ctmax ~ Species, TPC.pars))
##             Df Sum Sq Mean Sq F value Pr(>F)
## Species      3   98.2   32.73   1.772  0.177
## Residuals   26  480.2   18.47
summary(aov(e ~ Species, TPC.pars))
##             Df Sum Sq Mean Sq F value Pr(>F)
## Species      3 0.0021 0.00070    0.05  0.985
## Residuals   26 0.3612 0.01389
summary(aov(eh ~ Species, TPC.pars))
##             Df Sum Sq Mean Sq F value Pr(>F)
## Species      3  0.308  0.1028   0.296  0.828
## Residuals   17  5.910  0.3476               
## 9 observations deleted due to missingness
summary(aov(breadth ~ Species, TPC.pars))
##             Df Sum Sq Mean Sq F value Pr(>F)
## Species      3   9.92   3.306   1.888  0.156
## Residuals   26  45.52   1.751
summary(aov(topt ~ Species, TPC.pars))
##             Df Sum Sq Mean Sq F value Pr(>F)
## Species      3   3.11   1.037   0.623  0.606
## Residuals   26  43.26   1.664
d_fits <- mutate(d_fits, bootstrap = list(rep(NA, n())))

# run for loop to bootstrap each refitted model
for(i in 1:nrow(d_fits)){
  temp_data <- d_fits$data[[i]]
  temp_fit <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = temp_data,
                        start = coef(d_fits$sharpeschoolhigh[[i]]),
                        lower = get_lower_lims(temp_data$temp, temp_data$rate, model_name = 'sharpeschoolhigh_1981')-1,
                        upper = get_upper_lims(temp_data$temp, temp_data$rate, model_name = 'sharpeschoolhigh_1981')+1)
  boot <- Boot(temp_fit, method = 'residual')
  d_fits$bootstrap[[i]] <- boot
  rm(list = c('temp_fit', 'temp_data', 'boot'))
}
# get the raw values of each bootstrap
d_fits <- mutate(d_fits, output_boot = map(bootstrap, function(x) x$t))

# calculate predictions with a gnarly written function
d_fits <- mutate(d_fits, preds = map2(output_boot, data, function(x, y){
  temp <- as.data.frame(x) %>%
    drop_na() %>%
    mutate(iter = 1:n()) %>%
    group_by_all() %>%
    do(data.frame(temp = seq(min(y$temp), max(y$temp), length.out = 100))) %>%
    ungroup() %>%
    mutate(pred = sharpeschoolhigh_1981(temp, r_tref,e,eh,th, tref=28))
  return(temp)
}))

# select, unnest and calculate 95% CIs of predictions
boot_conf_preds <- select(d_fits, curve_id, preds) %>%
  unnest(preds) %>%
  group_by(curve_id, temp) %>%
  summarise(conf_lower = quantile(pred, 0.025),
            conf_upper = quantile(pred, 0.975),
            .groups = 'drop')

ggplot() +
  geom_line(aes(temp, .fitted), d_preds, col = 'blue') +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), boot_conf_preds, fill = 'blue', alpha = 0.3) +
  geom_point(aes(temp, rate), d, size = 2) +
  theme_bw(base_size = 12) +
  labs(x = 'Temperature (ºC)',
       y = 'Rate') +
  facet_wrap(~curve_id, ncol=4)

# GRANDIS

# grandis CURVE FIT
# P. grandis
d.grandis <- d %>% 
  filter(Species == "P. grandis")

#fit 
grandis.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = d,
                        iter = c(3,3,3,3),
                        start_lower = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') - 1,
                        start_upper = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') + 1,
                        lower = get_lower_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
                        upper = get_upper_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)

grandis.fit
## Nonlinear regression model
##   model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th,     tref = 28)
##    data: data
##  r_tref       e      eh      th 
##  0.1413  0.5775  3.4506 37.0533 
##  residual sum-of-squares: 0.1234
## 
## Number of iterations to convergence: 14 
## Achieved convergence tolerance: 1.49e-08
#generate the predicted data
grandis_new_data <- data.frame(temp = seq(min(d.grandis$temp), max(d.grandis$temp), 0.5))
grandis.preds <- augment(grandis.fit, newdata = grandis_new_data)

#calculate TPC parameters
grandis.TCP.res <- calc_params(grandis.fit) %>%
  mutate_all(round, 2)   # round 

grandis.TCP.res
##   rmax  topt ctmin ctmax    e   eh  q10 thermal_safety_margin thermal_tolerance
## 1 0.17 33.24 -5.65 47.29 0.51 1.11 1.93                 14.04             52.94
##   breadth skewness
## 1       9     -0.6
### Bootstrapping  curve fit    
# refit model using nlsLM
grandis.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = d.grandis,
                        start = coef(grandis.fit),
                        lower = get_lower_lims(d.grandis$temp, d.grandis$rate, model_name = 'sharpeschoolhigh_1981')-1,
                        upper = get_upper_lims(d.grandis$temp, d.grandis$rate, model_name = 'sharpeschoolhigh_1981')+1,
                        weights = rep(1, times = nrow(d.grandis)))

# bootstrap using case resampling
grandis.boot1 <- Boot(grandis.fit_nlsLM, method = 'case')

# look at the data
head(grandis.boot1$t)
##         r_tref         e       eh       th
## [1,] 0.1305408 0.4820667 3.920797 38.09578
## [2,] 0.1417562 0.5525771 3.148593 37.36422
## [3,] 0.1389140 0.5897576 3.334068 37.19541
## [4,] 0.1334456 0.4320001 4.887677 38.02246
## [5,] 0.1437546 0.6846489 3.129592 36.34727
## [6,] 0.1450414 0.5766447 2.596644 36.64574
# create predictions of each bootstrapped model
grandis.boot1_preds <- grandis.boot1$t %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(iter = 1:n()) %>%
  group_by_all() %>%
  do(data.frame(temp = seq(min(d.grandis$temp), max(d.grandis$temp), length.out = 100))) %>%
  ungroup() %>%
  mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))

# calculate bootstrapped confidence intervals
grandis.boot1_conf_preds <- group_by(grandis.boot1_preds, temp) %>%
  summarise(conf_lower = quantile(pred, 0.025),
            conf_upper = quantile(pred, 0.975)) %>%
  ungroup()

# plot bootstrapped CIs
grandis.CI.plot <- ggplot() +
  geom_line(aes(temp, .fitted), grandis.preds, col = 'darkgreen') +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), grandis.boot1_conf_preds, fill = 'darkgreen', alpha = 0.3) +
  geom_point(aes(temp, rate), d.grandis, size = 2, alpha = 0.5,col = 'darkgreen') +
  theme_bw(base_size = 12) +
  labs(x = 'Temperature (ºC)',
       y = 'Rate (µmol O2/cm2/s)')
grandis.CI.plot

P. meandrina

# meandrina CURVE FIT
# P. meandrina
d.meandrina <- d %>% 
  filter(Species == "P. meandrina")

#fit 
meandrina.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = d,
                        iter = c(3,3,3,3),
                        start_lower = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') - 1,
                        start_upper = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') + 1,
                        lower = get_lower_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
                        upper = get_upper_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)

meandrina.fit
## Nonlinear regression model
##   model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th,     tref = 28)
##    data: data
##  r_tref       e      eh      th 
##  0.1413  0.5775  3.4506 37.0533 
##  residual sum-of-squares: 0.1234
## 
## Number of iterations to convergence: 14 
## Achieved convergence tolerance: 1.49e-08
#generate the predicted data
meandrina_new_data <- data.frame(temp = seq(min(d.meandrina$temp), max(d.meandrina$temp), 0.5))
meandrina.preds <- augment(meandrina.fit, newdata = meandrina_new_data)

#calculate TPC parameters
meandrina.TCP.res <- calc_params(meandrina.fit) %>%
  mutate_all(round, 2)   # round 

meandrina.TCP.res
##   rmax  topt ctmin ctmax    e   eh  q10 thermal_safety_margin thermal_tolerance
## 1 0.17 33.24 -5.65 47.29 0.51 1.11 1.93                 14.04             52.94
##   breadth skewness
## 1       9     -0.6
### Bootstrapping  curve fit    
# refit model using nlsLM
meandrina.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = d.meandrina,
                        start = coef(meandrina.fit),
                        lower = get_lower_lims(d.meandrina$temp, d.meandrina$rate, model_name = 'sharpeschoolhigh_1981')-1,
                        upper = get_upper_lims(d.meandrina$temp, d.meandrina$rate, model_name = 'sharpeschoolhigh_1981')+1,
                        weights = rep(1, times = nrow(d.meandrina)))

# bootstrap using case resampling
meandrina.boot1 <- Boot(meandrina.fit_nlsLM, method = 'case')

# look at the data
head(meandrina.boot1$t)
##         r_tref         e       eh       th
## [1,] 0.1509585 0.6786125 3.527396 36.04762
## [2,] 0.1375984 0.5540254 3.779009 37.61593
## [3,] 0.1431038 0.6306717 3.714936 36.51225
## [4,] 0.1456467 0.5567787 3.469522 36.88690
## [5,] 0.1316263 0.3965938 5.982946 38.59573
## [6,] 0.1478011 0.5880170 3.168019 36.62537
# create predictions of each bootstrapped model
meandrina.boot1_preds <- meandrina.boot1$t %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(iter = 1:n()) %>%
  group_by_all() %>%
  do(data.frame(temp = seq(min(d.meandrina$temp), max(d.meandrina$temp), length.out = 100))) %>%
  ungroup() %>%
  mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))

# calculate bootstrapped confidence intervals
meandrina.boot1_conf_preds <- group_by(meandrina.boot1_preds, temp) %>%
  summarise(conf_lower = quantile(pred, 0.025),
            conf_upper = quantile(pred, 0.975)) %>%
  ungroup()

# plot bootstrapped CIs
meandrina.CI.plot <- ggplot() +
  geom_line(aes(temp, .fitted), meandrina.preds, col = 'orange') +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), meandrina.boot1_conf_preds, fill = 'orange', alpha = 0.3) +
  geom_point(aes(temp, rate), d.meandrina, size = 2, alpha = 0.5,col = 'orange') +
  theme_bw(base_size = 12) +
  labs(x = 'Temperature (ºC)',
       y = 'Rate (µmol O2/cm2/s)')
meandrina.CI.plot

P. tuahiniensis

# tuahiniensis CURVE FIT
# P. tuahiniensis
d.tuahiniensis <- d %>% 
  filter(Species == "P. tuahiniensis")

#fit 
tuahiniensis.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = d,
                        iter = c(3,3,3,3),
                        start_lower = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') - 1,
                        start_upper = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') + 1,
                        lower = get_lower_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
                        upper = get_upper_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)

tuahiniensis.fit
## Nonlinear regression model
##   model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th,     tref = 28)
##    data: data
##  r_tref       e      eh      th 
##  0.1413  0.5775  3.4506 37.0533 
##  residual sum-of-squares: 0.1234
## 
## Number of iterations to convergence: 14 
## Achieved convergence tolerance: 1.49e-08
#generate the predicted data
tuahiniensis_new_data <- data.frame(temp = seq(min(d.tuahiniensis$temp), max(d.tuahiniensis$temp), 0.5))
tuahiniensis.preds <- augment(tuahiniensis.fit, newdata = tuahiniensis_new_data)

#calculate TPC parameters
tuahiniensis.TCP.res <- calc_params(tuahiniensis.fit) %>%
  mutate_all(round, 2)   # round 

tuahiniensis.TCP.res
##   rmax  topt ctmin ctmax    e   eh  q10 thermal_safety_margin thermal_tolerance
## 1 0.17 33.24 -5.65 47.29 0.51 1.11 1.93                 14.04             52.94
##   breadth skewness
## 1       9     -0.6
### Bootstrapping  curve fit    
# refit model using nlsLM
tuahiniensis.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = d.tuahiniensis,
                        start = coef(tuahiniensis.fit),
                        lower = get_lower_lims(d.tuahiniensis$temp, d.tuahiniensis$rate, model_name = 'sharpeschoolhigh_1981')-1,
                        upper = get_upper_lims(d.tuahiniensis$temp, d.tuahiniensis$rate, model_name = 'sharpeschoolhigh_1981')+1,
                        weights = rep(1, times = nrow(d.tuahiniensis)))

# bootstrap using case resampling
tuahiniensis.boot1 <- Boot(tuahiniensis.fit_nlsLM, method = 'case')

# look at the data
head(tuahiniensis.boot1$t)
##         r_tref         e       eh       th
## [1,] 0.1301138 0.5511504 4.414331 37.26193
## [2,] 0.1420862 0.6200535 4.002197 36.88074
## [3,] 0.1225529 0.5349330 5.719269 38.16553
## [4,] 0.1357566 0.5963970 4.068060 37.37621
## [5,] 0.1295459 0.5457687 5.117105 37.60311
## [6,] 0.1365303 0.6017160 3.325320 36.68036
# create predictions of each bootstrapped model
tuahiniensis.boot1_preds <- tuahiniensis.boot1$t %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(iter = 1:n()) %>%
  group_by_all() %>%
  do(data.frame(temp = seq(min(d.tuahiniensis$temp), max(d.tuahiniensis$temp), length.out = 100))) %>%
  ungroup() %>%
  mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))

# calculate bootstrapped confidence intervals
tuahiniensis.boot1_conf_preds <- group_by(tuahiniensis.boot1_preds, temp) %>%
  summarise(conf_lower = quantile(pred, 0.025),
            conf_upper = quantile(pred, 0.975)) %>%
  ungroup()

# plot bootstrapped CIs
tuahiniensis.CI.plot <- ggplot() +
  geom_line(aes(temp, .fitted), tuahiniensis.preds, col = 'purple') +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), tuahiniensis.boot1_conf_preds, fill = 'purple', alpha = 0.3) +
  geom_point(aes(temp, rate), d.tuahiniensis, size = 2, alpha = 0.5,col = 'purple') +
  theme_bw(base_size = 12) +
  labs(x = 'Temperature (ºC)',
       y = 'Rate (µmol O2/cm2/s)')
tuahiniensis.CI.plot

P. verrucosa

# verrucosa CURVE FIT
# P. verrucosa
d.verrucosa <- d %>% 
  filter(Species == "P. verrucosa")

#fit 
verrucosa.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = d,
                        iter = c(3,3,3,3),
                        start_lower = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') - 1,
                        start_upper = get_start_vals(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981') + 1,
                        lower = get_lower_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
                        upper = get_upper_lims(d$temp, d$rate, model_name = 'sharpeschoolhigh_1981'),
                        supp_errors = 'Y',
                        convergence_count = FALSE)

verrucosa.fit
## Nonlinear regression model
##   model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th,     tref = 28)
##    data: data
##  r_tref       e      eh      th 
##  0.1413  0.5775  3.4506 37.0533 
##  residual sum-of-squares: 0.1234
## 
## Number of iterations to convergence: 14 
## Achieved convergence tolerance: 1.49e-08
#generate the predicted data
verrucosa_new_data <- data.frame(temp = seq(min(d.verrucosa$temp), max(d.verrucosa$temp), 0.5))
verrucosa.preds <- augment(verrucosa.fit, newdata = verrucosa_new_data)

#calculate TPC parameters
verrucosa.TCP.res <- calc_params(verrucosa.fit) %>%
  mutate_all(round, 2)   # round 

verrucosa.TCP.res
##   rmax  topt ctmin ctmax    e   eh  q10 thermal_safety_margin thermal_tolerance
## 1 0.17 33.24 -5.65 47.29 0.51 1.11 1.93                 14.04             52.94
##   breadth skewness
## 1       9     -0.6
### Bootstrapping  curve fit    
# refit model using nlsLM
verrucosa.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
                        data = d.verrucosa,
                        start = coef(verrucosa.fit),
                        lower = get_lower_lims(d.verrucosa$temp, d.verrucosa$rate, model_name = 'sharpeschoolhigh_1981')-1,
                        upper = get_upper_lims(d.verrucosa$temp, d.verrucosa$rate, model_name = 'sharpeschoolhigh_1981')+1,
                        weights = rep(1, times = nrow(d.verrucosa)))

# bootstrap using case resampling
verrucosa.boot1 <- Boot(verrucosa.fit_nlsLM, method = 'case')

# look at the data
head(verrucosa.boot1$t)
##         r_tref         e       eh       th
## [1,] 0.1525225 0.5886983 3.320467 36.52780
## [2,] 0.1511853 0.5745659 3.419156 36.81460
## [3,] 0.1487430 0.5655326 3.291840 37.15557
## [4,] 0.1507273 0.5694067 3.428851 37.05959
## [5,] 0.1489431 0.5968072 3.058033 36.80796
## [6,] 0.1438597 0.5573666 4.150847 37.09051
# create predictions of each bootstrapped model
verrucosa.boot1_preds <- verrucosa.boot1$t %>%
  as.data.frame() %>%
  drop_na() %>%
  mutate(iter = 1:n()) %>%
  group_by_all() %>%
  do(data.frame(temp = seq(min(d.verrucosa$temp), max(d.verrucosa$temp), length.out = 100))) %>%
  ungroup() %>%
  mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))

# calculate bootstrapped confidence intervals
verrucosa.boot1_conf_preds <- group_by(verrucosa.boot1_preds, temp) %>%
  summarise(conf_lower = quantile(pred, 0.025),
            conf_upper = quantile(pred, 0.975)) %>%
  ungroup()

# plot bootstrapped CIs
verrucosa.CI.plot <- ggplot() +
  geom_line(aes(temp, .fitted), verrucosa.preds, col = 'pink') +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), verrucosa.boot1_conf_preds, fill = 'pink', alpha = 0.3) +
  geom_point(aes(temp, rate), d.verrucosa, size = 2, alpha = 0.5,col = 'pink') +
  theme_bw(base_size = 12) +
  labs(x = 'Temperature (ºC)',
       y = 'Rate (µmol O2/cm2/s)')
verrucosa.CI.plot

 #set plot colors
cols <- c("grandis"="darkgreen","meandrina"="orange",  "tuahiniensis"="purple",  "verrucosa"="pink")
  
# plot data and model fit
TPC.plot <- ggplot(data=d, aes(x=temp)) +
   geom_point(aes(temp, rate, color="grandis"), d.grandis, size = 2, alpha = 0.5) +
   geom_point(aes(temp, rate, color="meandrina"), d.meandrina, size = 2, alpha = 0.5) +
   geom_point(aes(temp, rate, color="tuahiniensis"), d.tuahiniensis, size = 2, alpha = 0.5) +
   geom_point(aes(temp, rate, color="verrucosa"), d.verrucosa, size = 2, alpha = 0.5) +
   geom_line(aes(temp, .fitted), grandis.preds, col = 'darkgreen', size=2) +
   geom_line(aes(temp, .fitted), meandrina.preds, col = 'orange', size=2) +
   geom_line(aes(temp, .fitted), tuahiniensis.preds, col = "purple", size=2) +
   geom_line(aes(temp, .fitted), verrucosa.preds, col = "pink", size=2) +
   geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), grandis.boot1_conf_preds, fill = "darkgreen", alpha = 0.3) +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), meandrina.boot1_conf_preds, fill = "orange", alpha = 0.3) +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), tuahiniensis.boot1_conf_preds, fill = 'purple', alpha = 0.3) +
  geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), verrucosa.boot1_conf_preds, fill = 'pink', alpha = 0.3) +
   xlim(16,41)+
  #ylim(0,0.45)+
   scale_x_continuous(breaks=c(16,22,24,26,28,30,32,34,36,41))+
   theme_bw(base_size = 12) +
   scale_colour_manual(name="Species",values=cols)+
   theme(legend.position = "none",
         panel.border = element_blank(), panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
   labs(x = 'Temperature (ºC)',
        y = expression("Rate"~µmol~O[2] ~cm^{-2}~s^{-1}))

TPC.plot 

Bootstrap CI for all TPC parameters

broom::tidy(grandis.fit_nlsLM)
## # A tibble: 4 × 5
##   term   estimate std.error statistic  p.value
##   <chr>     <dbl>     <dbl>     <dbl>    <dbl>
## 1 r_tref    0.139   0.00609     22.8  1.16e-32
## 2 e         0.602   0.0742       8.12 1.82e-11
## 3 eh        3.18    0.464        6.86 3.07e- 9
## 4 th       37.0     0.726       51.0  3.91e-54
broom::tidy(meandrina.fit_nlsLM)
## # A tibble: 4 × 5
##   term   estimate std.error statistic  p.value
##   <chr>     <dbl>     <dbl>     <dbl>    <dbl>
## 1 r_tref    0.141   0.00573     24.6  4.91e-31
## 2 e         0.566   0.0753       7.51 6.19e-10
## 3 eh        3.77    0.632        5.97 1.93e- 7
## 4 th       37.1     0.621       59.7  5.39e-51
broom::tidy(tuahiniensis.fit_nlsLM)
## # A tibble: 4 × 5
##   term   estimate std.error statistic  p.value
##   <chr>     <dbl>     <dbl>     <dbl>    <dbl>
## 1 r_tref    0.132   0.00535     24.7  2.10e-34
## 2 e         0.562   0.0778       7.23 7.41e-10
## 3 eh        3.91    0.698        5.60 4.88e- 7
## 4 th       37.5     0.613       61.1  1.71e-58
broom::tidy(verrucosa.fit_nlsLM)
## # A tibble: 4 × 5
##   term   estimate std.error statistic  p.value
##   <chr>     <dbl>     <dbl>     <dbl>    <dbl>
## 1 r_tref    0.150   0.00576     26.1  1.44e-40
## 2 e         0.572   0.0621       9.21 3.76e-14
## 3 eh        3.30    0.447        7.38 1.41e-10
## 4 th       36.8     0.619       59.5  2.27e-67

#GRANDIS

#calculate all the TPC parameters
grandis.extra_params <- calc_params(grandis.fit_nlsLM) %>%
  pivot_longer(everything(), names_to =  'param', values_to = 'estimate')

#calculate CIs for all the TPC parameters
grandis.ci_extra_params <- Boot(grandis.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(grandis.fit_nlsLM)), R = 200, method = 'case') %>%
  confint(., method = 'perc') %>%
  as.data.frame() %>%
  rename(conf_lower = 1, conf_upper = 2) %>%
  rownames_to_column(., var = 'param') %>%
  mutate(method = 'case bootstrap')

#join the parameters and CIs  
grandis.ci_extra_params<- left_join(grandis.ci_extra_params, grandis.extra_params)
grandis.ci_extra_params$Treatment <- "P. grandis"

#meandrina

#calculate all the TPC parameters
meandrina.extra_params <- calc_params(meandrina.fit_nlsLM) %>%
  pivot_longer(everything(), names_to =  'param', values_to = 'estimate')

#calculate CIs for all the TPC parameters
meandrina.ci_extra_params <- Boot(meandrina.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(meandrina.fit_nlsLM)), R = 200, method = 'case') %>%
  confint(., method = 'perc') %>%
  as.data.frame() %>%
  rename(conf_lower = 1, conf_upper = 2) %>%
  rownames_to_column(., var = 'param') %>%
  mutate(method = 'case bootstrap')

#join the parameters and CIs  
meandrina.ci_extra_params<- left_join(meandrina.ci_extra_params, meandrina.extra_params)
meandrina.ci_extra_params$Treatment <- "P. meandrina"

#tuahiniensis

#calculate all the TPC parameters
tuahiniensis.extra_params <- calc_params(tuahiniensis.fit_nlsLM) %>%
  pivot_longer(everything(), names_to =  'param', values_to = 'estimate')

#calculate CIs for all the TPC parameters
tuahiniensis.ci_extra_params <- Boot(tuahiniensis.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(tuahiniensis.fit_nlsLM)), R = 200, method = 'case') %>%
  confint(., method = 'perc') %>%
  as.data.frame() %>%
  rename(conf_lower = 1, conf_upper = 2) %>%
  rownames_to_column(., var = 'param') %>%
  mutate(method = 'case bootstrap')

#join the parameters and CIs  
tuahiniensis.ci_extra_params<- left_join(tuahiniensis.ci_extra_params, tuahiniensis.extra_params)
tuahiniensis.ci_extra_params$Treatment <- "P. tuahiniensis"

#verrucosa

#calculate all the TPC parameters
verrucosa.extra_params <- calc_params(verrucosa.fit_nlsLM) %>%
  pivot_longer(everything(), names_to =  'param', values_to = 'estimate')

#calculate CIs for all the TPC parameters
verrucosa.ci_extra_params <- Boot(verrucosa.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(verrucosa.fit_nlsLM)), R = 200, method = 'case') %>%
  confint(., method = 'perc') %>%
  as.data.frame() %>%
  rename(conf_lower = 1, conf_upper = 2) %>%
  rownames_to_column(., var = 'param') %>%
  mutate(method = 'case bootstrap')

#join the parameters and CIs  
verrucosa.ci_extra_params<- left_join(verrucosa.ci_extra_params, verrucosa.extra_params)
verrucosa.ci_extra_params$Treatment <- "P. verrucosa"

#Join Species estimates and CIs

All_params <- rbind(grandis.ci_extra_params, meandrina.ci_extra_params,  tuahiniensis.ci_extra_params,  verrucosa.ci_extra_params)
All_params <- All_params %>% 
 mutate_if(is.numeric, round, 2)

#Plot all of the estimates
estimate.plots <- ggplot(All_params, aes(Treatment, estimate, color=Treatment)) +
  geom_point(size = 2) +
  scale_color_manual(name="Treatment", values=c("darkgreen","orange", "purple",  "pink"))+
  geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
  theme_bw() +
  facet_wrap(~param, scales = 'free_y') +
  scale_x_discrete('')

estimate.plots

# #filter to only the most relavent and well characterized parameters
 All_params <- All_params %>% 
   filter(!param=="ctmin") %>%
#   filter(!param=="ctmax") %>%
#   filter(!param=="eh") %>%
   filter(!param=="rmax") %>%
   filter(!param=="skewness") %>%
#   filter(!param=="topt") %>%
   filter(!param=="thermal_tolerance") #%>%
#   filter(!param=="q10") %>%
#   filter(!param=="e") %>%
#   filter(!param=="breadth")%>%
#   filter(!param=="thermal_safety_margin")

  
#view estimate plots
estimate.plots <- ggplot(All_params, aes(Treatment, estimate, color=Treatment)) +
  geom_point(size = 2) +
  scale_color_manual(name="Treatment", values=c("darkgreen","orange", "purple",  "pink"))+
  geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
  theme_bw() +
  labs(y = NULL)+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(), 
        strip.placement = "outside") +
  facet_wrap(~param, scales = 'free_y', nrow=2)+
             #labeller = as_labeller(c(e = "e (Energy)", eh = " eh (Energy)", rmax= "Rmax (~nmol~O[2] ~larva^{-1}~min^{-1})",topt="Topt (Temperature °C)")), strip.position = "left") +
  scale_x_discrete('')

estimate.plots

ggsave("output/TPC_estimates_SharpSchool_Respiration.pdf", estimate.plots, dpi=300, w=6, h=2, units="in")

#Plot Curve and Estimate Output

#generate a combined figure of TPCs and estimate plots
figure <- ggarrange(TPC.plot , estimate.plots,
                    labels = c("A", "B"),
                    ncol = 1, nrow = 2,
                    heights=c(1,0.5))
figure

ggsave("output/Respiration_TPC_and_estimates.pdf", figure, dpi=300, w=6, h=8, units="in")